  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"


begin 

 
  np     = 10 
  npa    = 4 
  nnx    = 2
  nny    = 2 
  
  fba = "BC_SP2_AVE_2004"             ;; Nov 
  fbb = "BC_SP2_CR-AVE_2006_WB57"     ;; Feb 
  fbc = "BC_SP2_TC4_2007_WB57"        ;; Aug  
  fbd = "BC_SP2_CARB_2008"            ;; June  
  fbe = "BC_SP2_ARCTAS-2008_DC8_04"   ;; Apr 
  fbf = "BC_SP2_ARCTAS-2008_P3B_04"   ;; Apr 
  fbg = "BC_SP2_ARCPAC-2008_WP-3D"    ;; Apr  
  fbh = "BC_SP2_ARCTAS-2008_DC8_JJ"   ;; Jun-Jul  
  fbi = "BC_SP2_ARCTAS-2008_P3B_JJ"   ;; Jun-Jul    
  fbj = "BC_SP2_ICE-L_2007"           ;; Nov-Dec  

  fbma = "HAM2R46H_2000_BC_profile" 
  fbmb = "E0001B_2000_BC_profile" 
  fbmc = "HAM2S30_2000_BC_profile" 

  fnb = "profile_"+fbb+".nc" 
  fnc = "profile_"+fbc+".nc" 
  fne = "profile_"+fbe+".nc" 
  fnf = "profile_"+fbf+".nc" 
  fng = "profile_"+fbg+".nc" 
  fnh = "profile_"+fbh+".nc" 
  fni = "profile_"+fbi+".nc" 
  fnj = "profile_"+fbj+".nc" 

  fma = "profile_BC_SP2_AVE-A_HGT" 
  fmb = "profile_BC_SP2_AVE-A_MEAN" 
  fmc = "profile_BC_SP2_AVE-B_HGT" 
  fmd = "profile_BC_SP2_AVE-B_MEAN" 
  fme = "profile_BC_SP2_CARB_H" 
  fmf = "profile_BC_SP2_CARB_med" 
  fmg = "profile_BC_SP2_CARB_std" 

  xaa = asciiread(fma,-1,"float") 
  yaa = asciiread(fmb,-1,"float") 
  xab = asciiread(fmc,-1,"float") 
  yab = asciiread(fmd,-1,"float") 
  xac = asciiread(fme,-1,"float") 
  yac = asciiread(fmf,-1,"float") 
  zac = asciiread(fmg,-1,"float") 
  xaa = xaa(::-1) 
  yaa = yaa(::-1) 
  xab = xab(::-1) 
  yab = yab(::-1) 
  yaa@_FillValue = -999 
  yab@_FillValue = -999 
  xaa@_FillValue = -999 
  xab@_FillValue = -999 

  xaa = where(ismissing(yaa),-999,xaa) 
  xab = where(ismissing(yab),-999,xab) 
 
  print("open file   : "+fnb) 
  print("open file   : "+fnc) 
  print("open file   : "+fne) 
  print("open file   : "+fnf) 
  print("open file   : "+fng) 
  print("open file   : "+fnh) 
  print("open file   : "+fni) 
  print("open file   : "+fnj) 

  fnma = fbma+".nc" 
  fnmb = fbmb+".nc" 
  fnmc = fbmc+".nc" 


  flb = addfile(fnb,"r") 
  flc = addfile(fnc,"r") 
  fle = addfile(fne,"r") 
  flf = addfile(fnf,"r") 
  flg = addfile(fng,"r") 
  flh = addfile(fnh,"r") 
  fli = addfile(fni,"r") 
  flj = addfile(fnj,"r")  
  flma = addfile(fnma,"r")  
  flmb = addfile(fnmb,"r")  
  flmc = addfile(fnmc,"r")  

  ;pfa = fla->profile 
  ;pfd = flb->profile 

  pfb = flb->profile 
  pfc = flc->profile 
  pfe = fle->profile 
  pff = flf->profile 
  pfg = flg->profile 
  pfh = flh->profile 
  pfi = fli->profile 
  pfj = flj->profile 
  BCPM = flma->BCPM
  BCPN = flmb->BCPM
  BCPO = flmc->BCPM

  BCPM(:,:,:) = BCPM(:,:,::-1) 
  BCPN(:,:,:) = BCPN(:,:,::-1) 
  BCPO(:,:,:) = BCPO(:,:,::-1) 

  nk = dimsizes(pfb(0,:)) 

  print("number of levels  :  "+nk ) 

;================================
; creat array for polygon 
;================================
printVarSummary(yac) 
printVarSummary(zac) 
  IP10 = 1 
  pga = new((/np,nk/),float)  
  pga(0,:) = pfb(IP10,:) 
  pga(0,:) = -999 ;;100
  pga(1,:) = pfb(IP10,:)  
  pga(2,:) = pfc(IP10,:)
  pga(3,:) = pfc(IP10,:)
  pga(3,:) = -999 ;;std 
  pga(4,:) = pfe(IP10,:)
  pga(5,:) = pff(IP10,:)
  pga(6,:) = pfg(IP10,:)
  pga(7,:) = pfh(IP10,:)
  pga(8,:) = pfi(IP10,:)
  pga(9,:) = pfj(IP10,:)

  IP90 = 2 
  pgb = new((/np,nk/),float)  
  pgb(0,:) = pfb(IP90,:)
  pgb(0,:) = -999 ;;100
  pgb(1,:) = pfb(IP90,:)
  pgb(2,:) = pfc(IP90,:)
  pgb(3,:) = -999 ;;std 
  pgb(4,:) = pfe(IP90,:)
  pgb(5,:) = pff(IP90,:)
  pgb(6,:) = pfg(IP90,:)
  pgb(7,:) = pfh(IP90,:)
  pgb(8,:) = pfi(IP90,:)
  pgb(9,:) = pfj(IP90,:)

  pgc = yac 
  pgc = yac - zac 
  pgd = yac 
  pgd = yac + zac 

;================================
; creat array for model A polygon
;================================

  IP10 = 1
  poa = new((/np,nk/),float)
  poa(0,:) = (/ BCPN(0,IP10,:) /)
  poa(1,:) = (/ BCPN(1,IP10,:) /)
  poa(2,:) = (/ BCPN(2,IP10,:) /)
  poa(3,:) = (/ BCPN(3,IP10,:) /)
  poa(4,:) = (/ BCPN(4,IP10,:) /)
  poa(5,:) = (/ BCPN(5,IP10,:) /)
  poa(6,:) = (/ BCPN(6,IP10,:) /)
  poa(7,:) = (/ BCPN(7,IP10,:) /)
  poa(8,:) = (/ BCPN(8,IP10,:) /)
  poa(9,:) = (/ BCPN(9,IP10,:) /)

  IP90 = 2
  pob = new((/np,nk/),float)
  pob(0,:) = (/ BCPN(0,IP90,:) /)
  pob(1,:) = (/ BCPN(1,IP90,:) /)
  pob(2,:) = (/ BCPN(2,IP90,:) /)
  pob(3,:) = (/ BCPN(3,IP90,:) /)
  pob(4,:) = (/ BCPN(4,IP90,:) /)
  pob(5,:) = (/ BCPN(5,IP90,:) /)
  pob(6,:) = (/ BCPN(6,IP90,:) /)
  pob(7,:) = (/ BCPN(7,IP90,:) /)
  pob(8,:) = (/ BCPN(8,IP90,:) /)
  pob(9,:) = (/ BCPN(9,IP90,:) /)

  ;poa(:,:) = where(ismissing(pga(:,:)),-999,poa) 
  ;pob(:,:) = where(ismissing(pgb(:,:)),-999,pob) 

;================================
; creat array for model B polygon
;================================

  IP10 = 1
  pha = new((/np,nk/),float)
  pha(0,:) = (/ BCPM(0,IP10,:) /)
  pha(1,:) = (/ BCPM(1,IP10,:) /)
  pha(2,:) = (/ BCPM(2,IP10,:) /)
  pha(3,:) = (/ BCPM(3,IP10,:) /)
  pha(4,:) = (/ BCPM(4,IP10,:) /)
  pha(5,:) = (/ BCPM(5,IP10,:) /)
  pha(6,:) = (/ BCPM(6,IP10,:) /)
  pha(7,:) = (/ BCPM(7,IP10,:) /)
  pha(8,:) = (/ BCPM(8,IP10,:) /)
  pha(9,:) = (/ BCPM(9,IP10,:) /)

  IP90 = 2
  phb = new((/np,nk/),float)
  phb(0,:) = (/ BCPM(0,IP90,:) /)
  phb(1,:) = (/ BCPM(1,IP90,:) /)
  phb(2,:) = (/ BCPM(2,IP90,:) /)
  phb(3,:) = (/ BCPM(3,IP90,:) /)
  phb(4,:) = (/ BCPM(4,IP90,:) /)
  phb(5,:) = (/ BCPM(5,IP90,:) /)
  phb(6,:) = (/ BCPM(6,IP90,:) /)
  phb(7,:) = (/ BCPM(7,IP90,:) /)
  phb(8,:) = (/ BCPM(8,IP90,:) /)
  phb(9,:) = (/ BCPM(9,IP90,:) /)

  ;pha(:,:) = where(ismissing(pga(:,:)),-999,pha)
  ;phb(:,:) = where(ismissing(pgb(:,:)),-999,phb)

;================================
; creat array for median obs 
;================================

  nln = 4 

  IP50 = 0
  pa = new((/np,nln,nk/),float)

  iln = 0 
  pa(0,iln,:) = -999
  pa(1,iln,:) = pfb(IP50,:)
  pa(2,iln,:) = pfc(IP50,:)
  pa(3,iln,:) = -999
  pa(4,iln,:) = pfe(IP50,:)
  pa(5,iln,:) = pff(IP50,:)
  pa(6,iln,:) = pfg(IP50,:)
  pa(7,iln,:) = pfh(IP50,:)
  pa(8,iln,:) = pfi(IP50,:)
  pa(9,iln,:) = pfj(IP50,:)


  iln = 1 
  pa(0,iln,:) = (/ BCPM(0,IP50,:) /)  
  pa(1,iln,:) = (/ BCPM(1,IP50,:) /)  
  pa(2,iln,:) = (/ BCPM(2,IP50,:) /) 
  pa(3,iln,:) = (/ BCPM(3,IP50,:) /) 
  pa(4,iln,:) = (/ BCPM(4,IP50,:) /) 
  pa(5,iln,:) = (/ BCPM(5,IP50,:) /) 
  pa(6,iln,:) = (/ BCPM(6,IP50,:) /) 
  pa(7,iln,:) = (/ BCPM(7,IP50,:) /) 
  pa(8,iln,:) = (/ BCPM(8,IP50,:) /) 
  pa(9,iln,:) = (/ BCPM(9,IP50,:) /)  

  ;pa(:,1,:) = where(ismissing(pa(:,0,:)),-999,pa(:,1,:))  

  iln = 2
  pa(0,iln,:) = (/ BCPN(0,IP50,:) /)
  pa(1,iln,:) = (/ BCPN(1,IP50,:) /)
  pa(2,iln,:) = (/ BCPN(2,IP50,:) /)
  pa(3,iln,:) = (/ BCPN(3,IP50,:) /)
  pa(4,iln,:) = (/ BCPN(4,IP50,:) /)
  pa(5,iln,:) = (/ BCPN(5,IP50,:) /)
  pa(6,iln,:) = (/ BCPN(6,IP50,:) /)
  pa(7,iln,:) = (/ BCPN(7,IP50,:) /)
  pa(8,iln,:) = (/ BCPN(8,IP50,:) /)
  pa(9,iln,:) = (/ BCPN(9,IP50,:) /)

  ;pa(:,2,:) = where(ismissing(pa(:,0,:)),-999,pa(:,2,:))

  iln = 3
  pa(0,iln,:) = (/ BCPO(0,IP50,:) /)
  pa(1,iln,:) = (/ BCPO(1,IP50,:) /)
  pa(2,iln,:) = (/ BCPO(2,IP50,:) /)
  pa(3,iln,:) = (/ BCPO(3,IP50,:) /)
  pa(4,iln,:) = (/ BCPO(4,IP50,:) /)
  pa(5,iln,:) = (/ BCPO(5,IP50,:) /)
  pa(6,iln,:) = (/ BCPO(6,IP50,:) /)
  pa(7,iln,:) = (/ BCPO(7,IP50,:) /)
  pa(8,iln,:) = (/ BCPO(8,IP50,:) /)
  pa(9,iln,:) = (/ BCPO(9,IP50,:) /)

  ;pa(:,3,:) = where(ismissing(pa(:,0,:)),-999,pa(:,3,:))



;================================
; 
;================================


  plotform = "eps" 
  plotform = "pdf" 
  plotname = "fig_BCPA" 


  C_tiXAxisString  = " BC mass (ng m~S~-3~N~)"  
  C_tiYAxisString  = " Altitude above sea level (km) "
  C_xyExplicitLegendLabels = (/" OBS"," HAM2"," HAM"," WET_MIX"/)

  C_pmLegendDisplayMode    = "Always"                           ; turn on legend
  C_pmLegendSide           = "Top"                              ; Change location of 
  L_lgPerimOn              = False
  R_lgPerimThicknessF      = 1.0 
  L_lgTitleOn              = False                              ; turn on legend title
  C_lgTitleString          = ""                                 ; create legend title
  R_lgTitleFontHeightF     = .009                               ; font of legend title
  R_lgLabelFontHeightF     = .016                               ; change font height
  R_pmLegendWidthF         = 0.12                               ; Change width and
  R_pmLegendHeightF        = 0.14                               ; height of legend.
  R_pmLegendParallelPosF   = 0.80                               ; move units right
  R_pmLegendOrthogonalPosF =-0.45                               ; move units down


  axisfontF = 0.020

  L_thin_border = 1

  markerS = 0.010
  markerF = 2.0


  
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
; define general wks
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  wks = gsn_open_wks(plotform,plotname)   


   gsn_merge_colormaps(wks,"StepSeq25","gsdtol")

  plot = new(np,graphic)  
  dumy = new(np,graphic)  
  dumx = new(np,graphic)  

;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
; main cycle of panels
;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


  ip = 0 


  ux = 1. * (/1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 /) 

  ux = ux -0.5 

  nhgt = dimsizes(ux) 

  umax = 1000. 
  umin = 1.e-2
  ymax = 18.
  ymin = 0.


  res = True
  res@gsnFrame = False 
  res@gsnDraw  = False 

  res@tmXTOn = False
  res@tmYROn = False

  res@tmYRBorderOn = False
  res@tmXTBorderOn = False
  
 

  res@tiXAxisString = C_tiXAxisString
  res@tiYAxisString = C_tiYAxisString
  res@gsnLeftString = "" 
  


  res@gsnLeftStringFontHeightF    = axisfontF
  res@gsnLeftStringOrthogonalPosF = 0.02
  res@gsnLeftStringParallelPosF   = 0.00
  res@gsnRightString = ""
 ;res@tiMainString   = ""   ; title





;==========================
; tune the shape and location of the plot
;==========================
 
; res@vpXF      = 0.15
; res@vpYF      = 0.9
  res@vpWidthF  = 0.4
  res@vpHeightF = 0.6

;==========================
; set max and min of axis vale
;==========================
 
  res@trXMinF  =  umin
  res@trXMaxF  =  umax
  res@trYMaxF  =  ymax
  res@trYMinF  =  ymin

  res@trXLog              = True
  res@trXLog              = True
  res@trXReverse          = False
  res@trYLog              = False
  res@trYReverse          = False

;==========================
; set axis font and axis-note font
;==========================
 
  if(L_thin_border.eq.0) then 

  res@tmBorderThicknessF     = 1.0

  res@tmYMajorGridThicknessF = 1.0
  res@tmXMajorGridThicknessF = 1.0

  res@tmXBMajorThicknessF    = 1.0
  res@tmXTMajorThicknessF    = 1.0
  res@tmYLMajorThicknessF    = 1.0
  res@tmYRMajorThicknessF    = 1.0
  res@tmXBMinorThicknessF    = 1.0
  res@tmXTMinorThicknessF    = 1.0
  res@tmYLMinorThicknessF    = 1.0
  res@tmYRMinorThicknessF    = 1.0

  end if


  res@tiYAxisFontHeightF     = axisfontF
  res@tmYLLabelFontHeightF   = axisfontF
  res@tmXBLabelFontHeightF   = axisfontF ; for using month name
  res@tiXAxisFontHeightF     = axisfontF
  res@tiYAxisFontHeightF     = axisfontF

  res@tmXTMajorLengthF = 0.008
  res@tmXBMajorLengthF = 0.008

  res@tmYLMajorLengthF = 0.012
  res@tmYRMajorLengthF = 0.012
  
  res@tmXTMajorOutwardLengthF = 0.008
  res@tmXBMajorOutwardLengthF = 0.008

  res@tmYLMajorOutwardLengthF = 0.012
  res@tmYRMajorOutwardLengthF = 0.012
  
 ;res@tmXBMode = "Explicit"                                       ; explicit labels
 ;res@tmXBValues        = (/1,10,100,1000,10000,50000/)                   ; location of labels 
 ;res@tmXBLabels        = (/1,10,100,1000,10000,50000/)                   ; labels themselves
  res@tmLabelAutoStride = True                                   ; nice stride on labels
  
  res@tmXBFormat = "#^se"                          ; nice stride on labels
  
  



  ;res@xyLineThicknesses = (/1.0,1.0,1.0/)
 


  do ip = 0, np-1  

  res@trXMinF  =  umin
  res@trXMaxF  =  umax
  res@trYMaxF  =  ymax
  res@trYMinF  =  ymin

  if(ip.eq.0) then 
     res@gsnLeftString  = "AVE Houston (Nov 2004) ~C~29~S~o~N~N-38~S~o~N~N, 88~S~o~N~W-98~S~o~N~W "
  end if 

  if(ip.eq.1) then
     res@gsnLeftString  = "CR-AVE (Feb 2006) ~C~1~S~o~N~S-10~S~o~N~N, 79~S~o~N~W-85~S~o~N~W"
     res@tiYAxisString = ""
  end if

  if(ip.eq.2) then
     res@gsnLeftString  = "TC4 (Aug 2007) ~C~2~S~o~N~N-12~S~o~N~N, 80~S~o~N~W-92~S~o~N~W"
  end if

  if(ip.eq.3) then
     res@gsnLeftString  = "CARB (Jun 2008) ~C~33~S~o~N~N-54~S~o~N~N, 105~S~o~N~W-127~S~o~N~W"
     res@tiYAxisString = ""
  end if


  ;C_gsnLeftString  = (/"AVE Houston","CR-AVE","TC4","CARB"/)  

  if(ip.eq.0) then 
     res@trXMinF  =  5.e-2 
  end if 
  if(ip.eq.3) then 
     res@trXMinF  =  1.e-1 
  end if 
 
;==========================
; set marker type
;==========================

  ;res@xyMonoMarkLineMode  = False
  res@xyMarkLineModes     = (/ "Lines", "Lines","Lines","Lines"/) ;;,"Markers","Markers","Markers"/) 
  ;res@xyMarkers           = (/ 6,6,6 /) 
  ;res@xyMarkerColors      = (/ "black","red","blue"/) 
  ;res@xyMarkerSizes       = (/ markerS, markerS, markerS, markerS, markerS, markerS /) 
  ;res@xyMarkerThicknesses = (/ 3.0, markerF, markerF,markerF, markerF, markerF /)  
  res@xyLineColors      = (/ "black","red","blue","green"/) 
  res@xyLineThicknesses   = (/ 2.0,2.0,2.0,2.0/)
  res@xyDashPatterns      = (/ 0,1,2,3/)
  print("ip = "+ip)

;==========================
; create respective legend
;==========================

  if(ip.le.10000) then 
  
  res@pmLegendDisplayMode    = "Always"                           ; turn on legend
  res@pmLegendSide           = "Top"                              ; change location of 
  res@pmLegendParallelPosF   = R_pmLegendParallelPosF             ; move units right
  res@pmLegendOrthogonalPosF = R_pmLegendOrthogonalPosF           ; move units down
  res@pmLegendWidthF         = R_pmLegendWidthF 
  res@pmLegendHeightF        = R_pmLegendHeightF 


  res@lgPerimOn              = L_lgPerimOn
  res@lgPerimThicknessF      = R_lgPerimThicknessF
  res@lgTitleOn              = L_lgTitleOn
  res@lgTitleString          = C_lgTitleString
  res@lgTitleFontHeightF     = R_lgTitleFontHeightF
  res@lgLabelFontHeightF     = R_lgLabelFontHeightF

  res@xyExplicitLegendLabels = C_xyExplicitLegendLabels
  res@lgRightMarginF         = 0.20
  res@lgLeftMarginF          = 0.20

  else

  res@pmLegendDisplayMode    = "NoCreate "                        ; turn off legend

  end if


  print("creat legend") 

;==========================
; create primitive plot
;==========================
  res@tfPolyDrawOrder    = "Predraw"               ; put line on top 

  plot(ip) = gsn_csm_xy(wks,pa(ip,:,:),ux,res)        ; draw y plot

  if(ip.eq.0) then  
    res@pmLegendDisplayMode    = "Never"                           ; turn on legend
;yaa = yaa / 2. + yab / 2. 
    res@gsnLeftString  = ""
    ploa = gsn_csm_xy(wks,yaa,xaa,res)   
    plob = gsn_csm_xy(wks,yab,xab,res)   
  end if 
 
  if(ip.eq.3) then 
    res@gsnLeftString  = ""
    res@pmLegendDisplayMode    = "Never"                           ; turn on legend
    ploc = gsn_csm_xy(wks,yac,xac,res)   
  end if 
 
  print("creat base plot") 

;==========================
; draw polygon 
;==========================

  xp    = new( (/2*nhgt/), float )
  yp    = new( (/2*nhgt/), float )

  do k=0,nhgt-1 
    yp(k)          = ux(k)
    xp(k)          = pga(ip,k)
    xp(2*nhgt-1-k) = pgb(ip,k) 
    yp(2*nhgt-1-k) = ux(k)
  end do

;;; i   = NhlNewColor(wks,1.0,0.71,0.76) ;fenhong
;;; j   = NhlNewColor(wks,0.64,0.71,0.8) ;huilan
;;;
;;; gsres                   = True  
;;; gsres@gsFillColor       = (/1.0,0.71,0.76/) 
;;; gsres@gsFillColor       = (/0.64,0.71,0.8/)

  i   = NhlNewColor(wks,1.00,0.71,0.76) ;fenhong
  j   = NhlNewColor(wks,0.64,0.71,0.80) ;huilan
  k   = NhlNewColor(wks,0.70,0.75,0.85) ;huilan

  gsres                   = True
 ;gsres@gsFillColor       = (/1.00,0.71,0.76/) 
 ;gsres@gsFillColor       = (/0.64,0.71,0.80/) 
 ;gsres@gsFillColor       = (/0.70,0.75,0.85/) 
  gsres@gsFillColor       = 54
 
  dumy(ip) = gsn_add_polygon (wks,plot(ip),xp,yp,gsres)


  delete(gsres) 

  print("creat polygon") 

;;  if(ip.eq.3) then 
;;
;;;==========================
;;; draw polygon
;;;==========================
;;
;;  nhga = 16 
;;  xpa    = new( (/2*nhga/), float )
;;  ypa    = new( (/2*nhga/), float )
;;
;;  do k=0,nhga-1
;;    ypa(k)          = xac(k)
;;    xpa(k)          = (/pgc(k) /) 
;;    xpa(2*nhga-1-k) = (/pgd(k)/) 
;;    ypa(2*nhga-1-k) = xac(k)
;;  end do
;;
;;  gsres                   = True
;;  gsres@gsFillColor       = 54
;;
;;  dumy(ip) = gsn_add_polygon (wks,plot(ip),xpa,ypa,gsres)
;;
;;  delete(gsres)
;;
;;  end if 

;==========================
; draw model polygon
;==========================

  xp    = new( (/2*nhgt/), float )
  yp    = new( (/2*nhgt/), float )

  do k=0,nhgt-1
    yp(k)          = ux(k)
    xp(k)          = (/ pha(ip,k)/) 
    xp(2*nhgt-1-k) = (/ phb(ip,k)/) 
    yp(2*nhgt-1-k) = ux(k)
  end do


  gsres                   = True
 ;gsres@gsFillColor       = (/1.00,0.71,0.76/)
 ;gsres@gsFillColor       = (/0.64,0.71,0.80/)
  gsres@gsFillColor       = (/0.70,0.75,0.85/)
  dumx(ip) = gsn_add_polygon (wks,plot(ip),xp,yp,gsres)

  delete(gsres) 

  print("creat polygon model ") 

  end do ;; loop ip

  overlay(plot(0),ploa) 
  overlay(plot(0),plob) 
  overlay(plot(3),ploc) 

;==========================
; create panel plot
;==========================
  resp = True
  resp@gsnPanelYWhiteSpacePercent = 5
  resp@gsnPanelXWhiteSpacePercent = 0

  gsn_panel(wks,plot(0:3),(/nnx,nny/),resp) 

  ;draw(plot)
  ;frame(wks)
   


 
end


